Analysis

BUSINESS QUESTIONS: Note: ignore the fact there is no question 1, this is due to a formatting error 2. a) Is there are a trend in injuries By Region, and is there different peaks of times of year per region? b) Analysis on the types of animals that are injured, this also by Region – is there a species that is more liable to injury in certain regions? - (cat/dog by region) c) What is the outcome? Does this differ by region?

1. Total call volume for complaint calls: How has this trended over time?
2. Is there a particular animal being called about the most?
3. Do particular suburbs have different type of complaint calls? Do they call about different animals? ((MAKE A LEAFLET MAP FOR THIS!))

4. Business Intelligence – using the insights you have found, can you predict how this might look for the upcoming year?

library(tidyverse)
library(tsibble)
library(forecast)
source("cleaning_script.R")

Rows: 31330 Columns: 7── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr (6): nature, animal_type, category, suburb, date_range, city
lgl (1): responsible_office
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Warning: Expected 2 pieces. Missing pieces filled with `NA` in 5643 rows [7, 11, 12, 18, 19, 25, 27, 29, 30, 31, 39, 44, 52, 54, 66, 67, 74, 75, 78, 81, ...].Rows: 42413 Columns: 5── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): Animal Type, Complaint Type, Date Received, Suburb, Electoral Division
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 664 Columns: 12── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): animal_type, outcome
dbl (10): year, ACT, NSW, NT, QLD, SA, TAS, VIC, WA, Total
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Note; summer in Australia seasons: Summer: December - February Autumn: March - May Winter: June - August Spring: September - November

also: The ‘wet season’ in Australia’s North: November - April

https://www.mdpi.com/2076-2615/8/7/100 will useful reading for later, talk about the need to reduce euthanasia or similar and the effect this has on people

Intro: Introduce the point of the talk, talk about Australian RSPCA, talk about how the data was gathered (using the information from the websites), the PURPOSE of this investigation (which is to help the RSPCA know which areas/animals to focus their efforts on), as I’m introducing the datasets I can introduce the two different cities. First, we can talk a short bit about Australia as a whole, it’s climate, the kind of animals etc. The point here is to really set the scene before diving too deep into facts and figures, as especially a non-technical audience this will help keep them engaged and make a more holistic presentation. In my opinion, it’s always good to zoom out and see the big picture, rather than getting lost in the myopia of some csv files.

Townsville Intro: Townsville is a city on the north-eastern coast of Queensland, Australia. With a population of 180,820 as of June 2018, it is the largest settlement in North Queensland; it is unofficially considered its capital. [note: put a map of Australia with Queensland and Townsville highlighted here. Talk a little bit about the population density, urbanisation, climate, types of animals that are common here etc. Show some photos of the area too]

Brisbane Intro: Exact same as above

3.

1. Total call volume for complaint calls: How has this trended over time?

First, let’s look at the Townsville animal complaints.

animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date_received, y = count)) +
  geom_line() +
  scale_x_date(date_breaks = "6 months", date_labels = "%b-%y")

This shows some seasonality and also an increase then a decline. Each summer (in December) the calls are much lower, rising again each Winter. Using geom_smooth, we get:

animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date_received, y = count)) +
  geom_smooth() +
  geom_point() +
  scale_x_date(date_breaks = "1 year", date_labels = "%y")

This shows a bit more clearly the general trend of call volume. From 2014 it steadily rises, peaking in 2017. Afterwards, it steadily declines to only slightly higher than where it started. It would be very difficult to say whether this trend will continue downward, go upward or stay relatively flat.

Note; if we have time or if it’s helpful, we will fix this graph so that quarters are properly displayed.

brisbane_complaints %>% 
  group_by(date) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date, y = count)) +
  geom_point() +
  geom_line() +
  scale_x_date(date_breaks = "3 months", date_labels = "%y")

Again, we see the seasonality of winter having more calls. The fact this is in both Brisbane and Townsville suggests a fairly general trend.

brisbane_complaints %>% 
  group_by(date) %>% 
  summarise(count = n()) %>% 
  ggplot(aes(x = date, y = count)) +
  geom_point() +
  geom_smooth() +
  scale_x_date(date_breaks = "1 year", date_labels = "20%y")

Other than declining a little over 2016, the number of calls sees a slow but steady increase towards 2020, being thousands more than it was in 2016 and 2017. So, the general trend is that the RSPCA are getting more complaints as time goes on. Now, this does not necessarily mean the line will continue to go up. We are also missing Q3 from 2016 so this skews the curve a little

3.b) Is there a particular animal being called about the most?

animal_complaints %>% 
  group_by(animal_type) %>% 
  ggplot(aes(x = animal_type)) +
  geom_bar()

The number of calls about dogs dwarf those about cats hugely. Let’s look at specific numbers:

animal_complaints %>% 
  group_by(animal_type) %>%
  count() %>% 
  summarise(n, 4094 / 38319)

So cats account for only 10% of the calls dogs do for Townsville! Let’s look at Brisbane:

brisbane_complaints %>% 
  group_by(type_of_animal) %>%
  count()

brisbane_complaints %>% 
  group_by(type_of_animal) %>%
  ggplot(aes(x = type_of_animal)) +
  geom_bar()

We have a lot more animal types here that we don’t in the previous data. We also have many calls about Attacks with no animal specified, many of which it is very likely they were involving dogs. However, we can’t say this for sure.

So that makes

4745 / 13334  # cats are 35% of the calls compared to dogs

[1] 0.3558572

# all other animals divided by dogs, leaving out unspecified (which we believe potentially contain a high proportion of dogs)
9457 / 13349

[1] 0.7084426

Most of the other animals have very small counts, interestingly foxes seem to make up a decent proportion of the calls regarding wild animals.

To answer the question however, it’s mainly dogs and cats, and especially dogs. Dogs are being called about more than any of the other animals combined (if we leave Unspecified to the side)

(Attack refers to the initial description of the complaint)

3. c) Do particular suburbs have different type of complaint calls? Do they call about different animals?

OK, so this one is difficult simply for the fact there is a huge amount of suburbs.

brisbane_complaints %>% 
  group_by(suburb) %>% 
  count()

# There is 192 suburbs! Certainly using a fill on a graph is not gonna work, neither is putting them on the x-axis on a bar graph.

animal_complaints %>% 
  group_by(suburb) %>% 
  count()

# 85 suburbs for Townsville

animal_complaints %>% 
  group_by(electoral_division) %>% 
  count()

# also 11 electoral divisions. Could we look at suburbs one electoral division at a time? Possibly

Idea: What about using leaflet to visualise types of complaint calls on a map?

This is more feasible with the Townsville datasets, as it has less suburbs and only 6 types of complaints. It may be necessary to try and wrangle the Brisbane data a little further to narrow down categories (for both type_of_animal and complaint_type)

FOR TOWNSVILLE WE HAVE: 85 different suburbs 6 complaint types 2 animal types

We want to break it down by suburb and complaint type, and then by suburb and animal type

animal_complaints %>% 
  ggplot(aes(x = suburb, fill = animal_type)) +
  geom_bar(position = "fill") +
  coord_flip() +
  scale_x_discrete(guide = guide_axis(n.dodge = 5))

animal_complaints %>% 
  ggplot(aes(x = suburb, fill =complaint_type)) +
  geom_bar(position = "fill") 

  #coord_flip() +
  #scale_x_discrete(guide = guide_axis(n.dodge = 5))

The challenge again, is having too many suburbs that we can’t get any useful information out of the graphs

animal_complaints %>% 
  group_by(suburb) %>% 
  count(sort = TRUE)

# let's drop any suburbs with less than 200 cases

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), animal_type) %>% 
  filter(count >= 500) %>% 
  ggplot(aes(x = suburb, y = count, fill = animal_type)) +
  geom_col() +
  coord_flip()

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

The highest count by quite a large margin is unallocated to a specific suburb.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), animal_type) %>% 
  filter(count >= 500) %>% 
  ggplot(aes(x = suburb, y = count, fill = animal_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

At a glance, there’s no big difference on animal types. All the suburbs have a large majority of dogs. Could we do a hypothesis test to see if there is a statistically significant difference?

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), animal_type) %>% 
  filter(count < 500 & count > 100) %>% 
  ggplot(aes(x = suburb, y = count, fill = animal_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

Looking at the lower end of the population, we see an outlier. Townsville City has way more cats than any of the other suburbs. In fact, it’s almost 50/50!

Now still looking at the Townsville dataset, we’ll break it down by complaint type:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500 & count <= 4000) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

Now for the lower end of the count:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 50) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

I’m not sure much can be gleamed with so many complaint_types. Let’s try focusing more specifically.

One thing to note here, is the vast variation in tolerance for noise. Most categories are consistent except this one. Cluden for example has a small percentage of noise complaints,while Bohle Plains has a huge percentage. This requires more investigation to figure out the root cause of this:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Noise" | (complaint_type == "Attack")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Noise" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Wandering" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count >= 500) %>% 
  filter(complaint_type == "Enclosure" | (complaint_type == "Noise")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

For the less than 500 greater than 50 groups:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>% 
  filter(complaint_type == "Noise" | (complaint_type == "Attack")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>%  
  filter(complaint_type == "Noise" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>% 
  filter(complaint_type == "Wandering" | (complaint_type == "Aggressive Animal")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() 

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type) %>% 
  filter(count < 500 & count > 100) %>% 
  #filter(complaint_type == "Enclosure" | (complaint_type == "Noise")) %>% 
  ggplot(aes(x = suburb, y = count, fill = complaint_type)) +
  geom_col(position = "fill") +
  coord_flip() +
  facet_wrap(~ suburb)

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type, animal_type) %>% 
  filter(count < 500 & count > 100) %>% 
  ggplot(aes(x = complaint_type, y = count, fill = complaint_type)) +
  geom_col() +
  facet_wrap(~ suburb)

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

From this graph, we can see that Hyde Park has a disproportionate amount of private impounds, while Bohle plains has more noise complaints.

Potential hypothesis test:

That the level of private impounds in Hyde Park being greater is statistically significant.

That the level of noise in Bohle Plains being greater is statistically significant.

animal_complaints %>% 
  group_by(suburb) %>% 
  count() %>% 
  filter(n < 500 & n > 100)

Let’s look at the suburbs with much higher complaints now:

animal_complaints %>% 
  group_by(suburb) %>% 
  summarise(count = n(), complaint_type, animal_type) %>% 
  filter(count >= 500 & count <4000) %>% 
  ggplot(aes(x = complaint_type, y = count, fill = complaint_type)) +
  geom_col() +
  facet_wrap(~ suburb)

`summarise()` has grouped output by 'suburb'. You can override using the `.groups` argument.

2. a) Are there any major differences between regions?

2. How do the animal outcomes differ over time?

3. Are there any noticeable anomalies in the data?

It’s impossible to say if there’s different peaks of times of yer per region, as we only have the data for the year as a whole. Unless we look at Suburb rather than region, and break this down by time. However, I think it’d be good to use the general Autralia (nationwide) data.

Proportions overall by region and by outcome. We can see that in general, regions have pretty similar proportions across all outcomes. This won’t be the case for every year though and trends over time, that’s the more interesting stuff.

This is the same thing as before, but for animal type. Same as before, not much discernible difference between the regions.

Next is just the sum ocurrences. We can see the regions that have way more animals injured/taken in by the RSPCA. This could be a population thing, maybe these regions just have more wildlife populations etc., or even a bigger RSPCA presence leading to more cases.

So, this is the data across time (from 1999-2018), with time on the x-axis and number of occurences on the y. Faceted into graphs based on outcome, with the coloured lines representing different regions.

What can we say about this? Well, QLD has had their transfers go way up since about 2009, and in general looks like an increase in cases for them. NSW had a big spike of Currently In Care around 2004. Euthanisation in general seems to be falling out of favour (which is actually true, the RSPCA tries to do it less for various reasons, one actually being the mental health of the staff and another being lots of euthanisations leads to a negative perception of animal shelters)

Exact same as before, but faceted on animal type rather than outcome.

QLD has had a lot more wildlife cases, steadily climbing since about 2006. NSW has had a big drop in cats recently. Dogs too. Horses have gone up every region it looks like. NT has had a huge spike in cats from 2015-2016.

This one has the same x and y, but is faceted based on region. The coloured lines now represent the outcomes.

WA in 2006 had a big increase in euthanizations. VIC has had theirs steadily decreasing over time, this is reflected also in NSW, TAS and ACT.

Same as before with faceting based on region, but now the coloured lines represent the type of animal the RSPCA is dealing with.

NT has had a big spike in dogs and cats in recent years. WA had a big spike in wildlife in 2005, corresponding with their spike in euthanizations - common sense and domain knowledge tell us these are likely highly correlated. VIC’s numbers in general seem to be decreasing over time.

Looking now specifically at euthanizations, faceted by region with the lines representing euthanizations.

Now we no longer facet by region.

OK so cats and dogs have way less euthanisations. Wildlife has gone way up! Horses are too small a sample to be meaningful. Other animals have spiked then decreased.

Looking at rehomes to euthanisations, we see cats have an upward trend with less euthanisations and more rehomes, with dogs, both have gone down. Perhaps with dogs there’s a confounding factor. In fact, it seems like the RSPCA is dealing with less dogs in general.

Euthanizations for wildlife have overall gone up. Queensland is skewing the data, they’ve euthanized way more! Is this necessarily a bad thing? It could indicate a bigger RSPCA presence in the area, and more vigilance in regards to controlling the wildlife.

General thoughts:

big change in 2006 for WA; cyclone Emma March 2006 could’ve released many wild animals leading to high euthanizations (this for 2.c) different territories have different approaches;

4. Business Intelligence – using the insights you have found, can you predict how this might look for the upcoming year?

animal_outcomes %>% 
  drop_na() %>% 
  group_by(region, outcome, year) %>%
  summarize(total_injuries = sum(number_of_occurences)) %>%
  ungroup() %>% 
  filter(outcome == "Euthanized") %>% 
  ggplot(aes(x = year, y = total_injuries, color = outcome)) +
  geom_line() +
  facet_wrap(~region, scales = "free_y") +
  labs(x = "Date", y = "Number of injuries", color = "Outcome")

`summarise()` has grouped output by 'region', 'outcome'. You can override using the `.groups` argument.

Next attempt for a model

I want to predict the total cases for Australia. We won’t break it down by region. So first, let’s get our summary statistics and put it in a new variable:

animal_outcomes %>% 
  drop_na() %>% 
  group_by(outcome, animal_type, year, region) %>%
  summarize(total_injuries = sum(number_of_occurences)) %>%
  ungroup() %>% 
  filter(animal_type == "Cats" & outcome == "Euthanized") %>% 
  ggplot(aes(x = year, y = total_injuries, color = outcome)) +
  geom_line() +
  facet_wrap(~region, scales = "free_y") +
  labs(x = "Date", y = "Number of injuries", color = "Outcome", title = "Cats")

`summarise()` has grouped output by 'outcome', 'animal_type', 'year'. You can override using the `.groups` argument.

animal_outcomes %>% 
  drop_na() %>% 
  group_by(outcome, animal_type, year) %>%
  summarize(total_injuries = sum(number_of_occurences)) %>%
  ungroup() %>% 
  filter(animal_type == "Cats" & outcome == "Euthanized") %>% 
  ggplot(aes(x = year, y = total_injuries, color = outcome)) +
  geom_line() +
  #facet_wrap(~region, scales = "free_y") +
  labs(x = "Date", y = "Number of injuries", color = "Outcome", title = "Cats")

`summarise()` has grouped output by 'outcome', 'animal_type'. You can override using the `.groups` argument.

wildlife_tsbl <- animal_outcomes %>% 
  drop_na() %>% 
  filter(animal_type == "Cats" & outcome == "Euthanized") %>% 
  group_by(year) %>% 
  summarise(total_euth = sum(number_of_occurences)) %>% 
  ungroup() %>% 
  as_tsibble(key = year, index = total_euth)

# Fit ARIMA model
wildlife_arima <- wildlife_tsbl %>% 
  model(total_euth = ARIMA(total_euth))

Warning: 20 errors (1 unique) encountered for total_euth
[20] attempt to select less than one element in get1index

# Generate forecasts
wildlife_forecast <- wildlife_tsbl %>%
  model(ARIMA(total_euth)) %>% 
  forecast(h = 36)

Error in FUN(X[[i]], ...) : object 'total_euth' not found

library(tidyverse)
library(fable)
library(tsibble)
library(tsibbledata)
# We also need to use the "urca" package to be able to use the ARIMA model
library(urca)

library(tsibble)
library(fable)

wildlife_tsbl <- animal_outcomes %>% 
  drop_na() %>% 
  filter(animal_type == "Cats" & outcome == "Euthanized") %>% 
  group_by(year) %>% 
  summarise(total_euth = sum(number_of_occurences)) %>% 
  ungroup() %>%
  mutate(date = year %>% as.Date(paste0("-01-01"))) %>% 
  as_tsibble(key = year, index = date)

wildlife_arima <- wildlife_tsbl %>% 
  model(ARIMA(total_euth))

wildlife_forecast <- wildlife_arima %>% 
  forecast(h = 36)

wildlife_forecast

NA
NA

wildlife_forecast %>% 
  autoplot() + 
  scale_y_continuous(limits = c(0, max(wildlife_tsbl$total_euth, na.rm = TRUE) * 1.2))

library(fable)
library(tsibbledata)
library(dplyr)
library(ggplot2)


# Create tsibble
wildlife_tsbl <- animal_outcomes %>% 
  drop_na() %>% 
  filter(animal_type == "Cats" & outcome == "Euthanized") %>% 
  group_by(year) %>% 
  summarise(total_euth = sum(number_of_occurences)) %>% 
  ungroup() %>% 
  as_tsibble(index = year)

# Fit ARIMA model
wildlife_arima <- wildlife_tsbl %>% 
  model(ARIMA(total_euth))

Warning: 1 error encountered for ARIMA(total_euth)
[1] .data contains implicit gaps in time. You should check your data and convert implicit gaps into explicit missing values using `tsibble::fill_gaps()` if required.

# Generate forecasts
wildlife_fcst <- wildlife_arima %>% 
  forecast(h = 36)

# Add source column to forecast
wildlife_fcst %>%
  mutate(source = "forecast") %>%
  select(source, year, .mean) %>%
  bind_rows(wildlife_tsbl %>% as_tibble() %>% mutate(source = "actual")) %>%
  
  # Plot the actuals and forecasts
  ggplot(aes(x = year, y = .mean, color = source)) +
  geom_line() +
  labs(y = "Total Euthanized Cats") +
  scale_color_manual(values = c("actual" = "black", "forecast" = "red")) +
  theme_minimal()

Error in `build_fable()`:
! Can't subset columns that don't exist.
✖ Column `total_euth` doesn't exist.
Backtrace:
 1. ... %>% ggplot(aes(x = year, y = .mean, color = source))
 6. fabletools:::select.fbl_ts(., source, year, .mean)
 7. fabletools:::restore_fable(res, .data)
 8. fabletools:::build_fable(data, response = response_vars(template), distribution = !!distribution_var(template))

wildlife_ts

Time Series:
Start = 1 
End = 20 
Frequency = 1 
 [1] 31202 32901 31009 30449 35934 31941 37010 33343 42731 39495 38642 37177 24651 19464
[15] 15491 17398 16205 14563 12610 11740

library(tsibble)
library(fable)

# Convert ts object to tsibble
wildlife_tsbl <- as_tsibble(wildlife_ts, index = "year")

# Fit ARIMA model to tsibble
wildlife_arima <- wildlife_tsbl %>% 
  model(ARIMA(value))

# Generate forecasts
wildlife_forecast <- wildlife_arima %>% 
  forecast(h = 36)

# Convert forecast to a tsibble
wildlife_fcst <- as_tsibble(wildlife_forecast, key = "year")

# Plot the forecast
wildlife_fcst %>%
  ggplot(aes(x = year, y = .mean, color = source)) +
  geom_line()

Error in `geom_line()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 1st layer.
Caused by error in `compute_aesthetics()`:
! Aesthetics are not valid data columns.
✖ The following aesthetics are invalid:
✖ `x = NULL`
✖ `colour = NULL`
ℹ Did you mistype the name of a data column or forget to add `after_stat()`?
Backtrace:
  1. base (local) `<fn>`(x)
  2. ggplot2:::print.ggplot(x)
  4. ggplot2:::ggplot_build.ggplot(x)
  5. ggplot2:::by_layer(...)
 12. ggplot2 (local) f(l = layers[[i]], d = data[[i]])
 13. l$compute_aesthetics(d, plot)
 14. ggplot2 (local) compute_aesthetics(..., self = self)

library(tsibble)
library(fable)
library(ggplot2)

# Convert ts object to tsibble
wildlife_tsbl <- as_tsibble(wildlife_ts, index = "year")

# Fit ARIMA model to tsibble
wildlife_arima <- wildlife_tsbl %>% 
  model(ARIMA(value))

# Generate forecasts
wildlife_forecast <- wildlife_arima %>% 
  forecast(h = 36)

# Convert forecast to a tsibble
wildlife_fcst <- as_tsibble(wildlife_forecast, key = "year")

# Plot the forecast
wildlife_fcst %>%
  ggplot(aes(x = index, y = .mean, color = .model)) +
  geom_line()

NA
NA

names(wildlife_fcst)

[1] ".model" "index"  "value"  ".mean" 

# grouping the data and creating a variable
animal_complaints_grouped <- animal_complaints %>% 
  group_by(date_received) %>% 
  summarise(count = n()) %>% 
  arrange(date_received)

myts <- ts(animal_complaints_grouped$count, frequency = 12, start = c(year(animal_complaints_grouped$date_received[1]), month(animal_complaints_grouped$date_received[1])))

myforecast <- forecast(auto.arima(myts), h = 24)

Error in switch(names(interval), day = "days", hour = "hours", minute = "mins",  : 
  EXPR must be a length 1 vector

myts_forecast <- function(myts, h, interval) {
  library(forecast)
  
  # Set the appropriate time interval for the forecast horizon
  interval <- switch(interval,
                     "day" = "days",
                     "hour" = "hours",
                     "minute" = "mins",
                     "second" = "secs")
  
  # Create the forecast
  myforecast <- forecast(myts, h = 24, level = c(80, 95), 
                       xreg = mydummies, 
                       lambda = NULL, biasadj = TRUE,
                       fan = FALSE, 
                       interval = "week")

  
  # Print the forecast table
  print(myforecast)
  
  # Plot the forecast
  autoplot(myforecast)
}

print(interval)

function (start = NULL, end = NULL, tzone = tz(start)) 
{
    if (is.character(start) && is.null(end)) {
        return(parse_interval(start, tzone))
    }
    if (length(start) == 0 || length(end) == 0) {
        if (is.null(start) && missing(tzone)) {
            tzone <- "UTC"
        }
        start <- POSIXct(tz = tzone)
        return(new("Interval", numeric(), start = start, tzone = tzone))
    }
    if (is.Date(start)) 
        start <- date_to_posix(start)
    if (is.Date(end)) 
        end <- date_to_posix(end)
    force(tzone)
    start <- as_POSIXct(start, tzone)
    end <- as_POSIXct(end, tzone)
    span <- as.numeric(end) - as.numeric(start)
    starts <- start + rep(0, length(span))
    if (tzone != tz(starts)) 
        starts <- with_tz(starts, tzone)
    new("Interval", span, start = starts, tzone = tzone)
}
<bytecode: 0x7fa7b16da468>
<environment: namespace:lubridate>

print(class(interval))

[1] "function"

print(length(interval))

[1] 1

# Create a sequence of dates for the next 24 months
dates <- seq(as.Date("2023-03-01"), by = "month", length.out = 24)

# Create a time series object with random values for demonstration purposes
values <- rnorm(24)
myts <- ts(values, start = c(year(dates[1]), month(dates[1])), frequency = 12)

# Plot the time series
plot(myts)

# Forecast the next 24 months
library(forecast)
myforecast <- forecast(myts, h = 24)

Error in UseMethod("forecast") : 
  no applicable method for 'forecast' applied to an object of class "ts"

Another, FINAL attempt

This is quite difficult indeed. Let’s try go slowly using the notes, build it bit by bit.

So all we wanna do, is plot the next few years, and see what it will look like on the graph. Nothing needs to be too fancy. We will keep it simple, create a simple model. Then plot this. OK.

So first, it needs to be tsibble.

It has now been converted to a tsibble, with date_received as the index. As it contains only a single time series, a key should not be necessary.

library(tidyverse)
library(fable)
library(tsibble)
library(tsibbledata)
# We also need to use the "urca" package to be able to use the ARIMA model
library(urca)

We load in all packages just to be 100% sure they are loaded.

# next, we will use autoplot
autoplot(animal_complaints_tsibble)

Plot variable not specified, automatically selected `.vars = count`

Our data is there. Now, we will fit a model.

library(purrr)

fit <- animal_complaints_tsibble %>%
  model(
    snaive = SNAIVE(count),
    mean_model = MEAN(count),
    arima = ARIMA(count)
  )
fit

# fit %>% 
#   dmap(unlist)

Here is our mable, also called model table. Each cell corresponds to a fitted model. Now, we will calculate some forecasts using our model.

forecast_1 <- fit %>%
  fabletools::forecast(h = 36)
forecast_1

OK, looking good so far.

forecast_1 %>%
  autoplot(animal_complaints_tsibble)

Let’s make our model look pretty and have less clutter

forecast_1 %>%
  filter(.model == "arima") %>% 
  autoplot(animal_complaints_tsibble, level = 95) +
  scale_x_yearmonth(date_breaks = "6 months", date_labels = "%b %y") + 
  labs(title = "Townsville Complaint Calls with Prediction",
       x = "Time",
       y = "Number of Calls") +
  theme_economist() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, family = "mono", face = "bold", size = 7),
        #axis.ticks.x = 
        axis.title.y = element_text(vjust = 3.2),
        title = element_text(family = "mono", face = "bold"),
        plot.title = element_text(hjust = 0.4),
        legend.title = element_blank(),
        legend.position = "none") +
        annotate("rect", xmin = as.Date("2014-06-01"), xmax = as.Date("2014-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2015-06-01"), xmax = as.Date("2015-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2016-06-01"), xmax = as.Date("2016-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2017-06-01"), xmax = as.Date("2017-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2018-06-01"), xmax = as.Date("2018-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2019-06-01"), xmax = as.Date("2019-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2020-06-01"), xmax = as.Date("2020-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2021-06-01"), xmax = as.Date("2021-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2022-06-01"), xmax = as.Date("2022-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3) +
        annotate("rect", xmin = as.Date("2023-06-01"), xmax = as.Date("2023-08-01"), ymin = -Inf, ymax = Inf, fill = "gray70", alpha = 0.3)

NA

So I was also tasked with creating a model, to predict how things might look for the next few years. I focused specifically on the Townsville data, as this was most suitable for forecasting. I used something called the “ARIMA” model, the inner workings of which are quite complicated, but basically it predicts the future based on the past, and is quite good when it comes to seasonal data like I have here. The blue line is the prediction, and the light blue area you see highlighted around it is the confidence interval. Think of this like a margin of error, and it basically means that 95% of the time, our model will predict results in this area.

This model suggests that roughly the same pattern will continue as before, with fairly steady volume of calls and high peaks in Winter. But bear in mind this model is quite limited and shouldn’t be taken as fact; for example, it doesn’t take into account natural disasters and other events which could cause an increase in complaint calls.

Documentation note

Fix or discuss the error of NT and NSW swapping for 2016-2018